//#include "../DWMRICore/Point3F.h"
//#include "../DWMRICore/Point3I.h"
//#include "../DWMRICore/VectorField.h"
//#include "../DWMRICore/Fiber.h"
//#include "../DWMRICore/Scalar.h"
//
//#include <stdio.h>
//#include <stdlib.h>
//#include <string.h>
//
//#include <vector>
//
//using namespace std;
//
//char *src_pathname = NULL;
//char *edge_pathname = NULL;
//char *dst_pathname = NULL;
//
//CVectorField *vec = NULL;
//CScalar *edge = NULL;
//int w;
//int h;
//int d;
//
//const float threshold1 = 0.9993f;
//const float threshold2 = 0.866f;
//const float edge_theshold = 200.0f;
//
//struct Element{
//	CVector3F dir;
//	CPoint3I pos;
//	float weight;
//};
//
//void PrintUsage()
//{
//	printf("Usage: RegionGrowing -src input_vector_file -edge input_edge_file -dst output_vector_file\n");
//}
//
//int ParseArguments(const int argc, char *argv[])
//{
//	for (int i = 1; i < argc; ++i) {
//		if (strcmp(argv[i], "-src") == 0) {
//			src_pathname = argv[++i];
//		} else if (strcmp(argv[i], "-edge") == 0) {
//			edge_pathname = argv[++i];
//		} else if (strcmp(argv[i], "-dst") == 0) {
//			dst_pathname = argv[++i];
//		} else {
//			return -1;
//		}
//	}
//
//	if (src_pathname == NULL || dst_pathname == NULL || edge_pathname == NULL)
//		return -2;
//
//	return 0;
//}
//
//void ReadVectorFields(vector<vector<CVector3F> > &vectorFields, int &w, int &h, int &d)
//{
//	FILE *fp = fopen(src_pathname, "rb");
//
//	fread(&w, sizeof(int), 1, fp);
//	fread(&h, sizeof(int), 1, fp);
//	fread(&d, sizeof(int), 1, fp);
//
//	vectorFields.resize(w*h*d);
//
//	for (int i = 0; i < w * h * d; ++i) {
//		int count;
//		fread(&count, sizeof(int), 1, fp);
//
//		for (int j = 0; j < count; ++j) {
//			float v[3];
//			fread(v, sizeof(float), 3, fp);
//			vectorFields.at(i).push_back(CVector3F(v[0], v[1], v[2]));
//		}
//	}
//
//	fclose(fp);
//}
//
////void Inseart(Element &s, CVector3F &dir, vector<Element> &maxHeap, CPoint3I &pos)
////{
////	float weight = fabs(InnerProduct(s.dir, dir));
////
////	// inseart it to the max heap
////	uint p;
////	uint left = 0;
////	uint right = maxHeap.size();
////
////	while ((right - left) > 10) {
////		int middle = (left + right) / 2;
////		if (maxHeap.at(middle).weight > weight)
////			right = middle;
////		else
////			left = middle;
////	}
////
////	for (p = left; p < right; ++p) {
////		if (maxHeap.at(p).weight > weight)
////			break;
////	}
////	
////	Element e;
////	e.dir = dir;
////	e.pos = pos;
////	e.weight = weight;
////
////	maxHeap.insert(maxHeap.begin()+p, e);
////}
//
//void Inseart(Element &s, CVector3F &dir, vector<Element> &maxHeap, CPoint3I &pos)
//{
//	Element e;
//	e.dir = dir;
//	e.pos = pos;
//	maxHeap.push_back(e);
//}
//
//void RegionGrowing(CPoint3I &seed, float *res)
//{
//	memset(res, 0, sizeof(float)*w*h*d);
//
//	vector<Element> maxHeap;
//	int *flags = new int[w*h*d];
//	memset(flags, 0, sizeof(int)*w*h*d);
//
//	CVector3F seed_dir = vec->GetVectorFast(seed.x, seed.y, seed.z);
//
//	Element e;
//	e.dir = vec->GetVectorFast(seed.x, seed.y, seed.z);
//	e.pos = seed;
//	e.weight = 1.0f;
//	maxHeap.push_back(e);
//	
//	int index = (seed.z * h + seed.y) * w + seed.x;
//	flags[index] = 1;
//
//	int done = 1;
//	int last_done = 1;
//	while(maxHeap.size() > 0) {
//		Element s = maxHeap.back();
//		maxHeap.pop_back();
//
//		index = (s.pos.z * h + s.pos.y) * w + s.pos.x;
//		CPoint3I pos;
//		pos = s.pos;
//
//		if (pos.z < 140 || pos.z > 390)
//			continue;
//
//		res[index] =  1.0f;
//
//		CVector3F temp_dir;
//		float edge_strength;
//
//		temp_dir = vec->GetVectorFast(pos.x-1, pos.y, pos.z);
//		edge_strength = edge->GetScalarFast(pos.x-1, pos.y, pos.z);
//		if (pos.x > 0 && flags[index-1] == 0 && fabs(InnerProduct(temp_dir, s.dir)) > threshold1 && 
//			fabs(InnerProduct(temp_dir, seed_dir)) > threshold2 && edge_strength < edge_theshold) {
//			Inseart(s, temp_dir, maxHeap, CPoint3I(pos.x-1, pos.y, pos.z));
//			flags[index-1] = 1;
//			done++;
//		}
//
//		temp_dir = vec->GetVectorFast(pos.x, pos.y-1, pos.z);
//		edge_strength = edge->GetScalarFast(pos.x, pos.y-1, pos.z);
//		if (pos.y > 0 && flags[index-w] == 0  && fabs(InnerProduct(temp_dir, s.dir)) > threshold1 && 
//			fabs(InnerProduct(temp_dir, seed_dir)) > threshold2 && edge_strength < edge_theshold) {
//			Inseart(s, temp_dir, maxHeap, CPoint3I(pos.x, pos.y-1, pos.z));
//			flags[index-w] = 1;
//			done++;
//		}
//
//		temp_dir = vec->GetVectorFast(pos.x, pos.y, pos.z-1);
//		edge_strength = edge->GetScalarFast(pos.x, pos.y, pos.z-1);
//		if (pos.z > 0 && flags[index-w*h] == 0 && fabs(InnerProduct(temp_dir, s.dir)) > threshold1 && 
//			fabs(InnerProduct(temp_dir, seed_dir)) > threshold2 && edge_strength < edge_theshold) {
//			Inseart(s, temp_dir, maxHeap, CPoint3I(pos.x, pos.y, pos.z-1));
//			flags[index-w*h] = 1;
//			done++;
//		}
//
//		temp_dir = vec->GetVectorFast(pos.x+1, pos.y, pos.z);
//		edge_strength = edge->GetScalarFast(pos.x+1, pos.y, pos.z);
//		if (pos.x < w-1 && flags[index+1] == 0 && fabs(InnerProduct(temp_dir, s.dir)) > threshold1 && 
//			fabs(InnerProduct(temp_dir, seed_dir)) > threshold2 && edge_strength < edge_theshold) {
//			Inseart(s, temp_dir, maxHeap, CPoint3I(pos.x+1, pos.y, pos.z));
//			flags[index+1] = 1;
//			done++;
//		}
//
//		temp_dir = vec->GetVectorFast(pos.x, pos.y+1, pos.z);
//		edge_strength = edge->GetScalarFast(pos.x, pos.y+1, pos.z);
//		if (pos.y < h-1 && flags[index+w] == 0 && fabs(InnerProduct(temp_dir, s.dir)) > threshold1 && 
//			fabs(InnerProduct(temp_dir, seed_dir)) > threshold2 && edge_strength < edge_theshold) {
//			Inseart(s, temp_dir, maxHeap, CPoint3I(pos.x, pos.y+1, pos.z));
//			flags[index+w] = 1;
//			done++;
//		}
//
//		temp_dir = vec->GetVectorFast(pos.x, pos.y, pos.z+1);
//		edge_strength = edge->GetScalarFast(pos.x, pos.y, pos.z+1);
//		if (pos.z < d-1 && flags[index+w*h] == 0 && fabs(InnerProduct(temp_dir, s.dir)) > threshold1 && 
//			fabs(InnerProduct(temp_dir, seed_dir)) > threshold2 && edge_strength < edge_theshold) {
//			Inseart(s, temp_dir, maxHeap, CPoint3I(pos.x, pos.y, pos.z+1));
//			flags[index+w*h] = 1;
//			done++;
//		}
//
//		if ((done - last_done) > 1000) {
//			float ff = (float)done / (float)(w*h*d) * 100.0f;
//			printf("\r%.3f%%", ff);
//			last_done = done;
//		}
//	}
//
//	delete[] flags;
//}
//
//int main(int argc, char *argv[])
//{
//	if (ParseArguments(argc, argv) != 0) {
//		PrintUsage();
//		return 0;
//	}
//
//	vec = new CVectorField();
//	vec->OpenVectorFieldFile(src_pathname);
//	
//	edge = new CScalar();
//	edge->OpenScalarFile(edge_pathname);
//
//	w = vec->GetWidth();
//	h = vec->GetHeight();
//	d = vec->GetDepth();
//
//	CPoint3I seed(231, 241, 295);
//
//
//	float *result = new float[w*h*d];
//
//	RegionGrowing(seed, result);
//
//
//	CScalar *res = new CScalar();
//	res->CreateScalar(w, h, d, result);
//	res->SaveScalarFile(dst_pathname);
//
//	SafeDelete(res);
//	SafeDelete(vec);
//	SafeDelete(edge);
//
//	return 0;
//}
//
////int main()
////{
////	char *fiber_pathname = "\\\\nabla.cs.purdue.edu\\visdata\\Ziang\\CoherenceStructure\\Datasets\\HumanConnectome\\792564_processed\\NEW_DIT_Fibers\\laptop_fibers.txt";
////	
////	char **confPathname = NULL;
////	char **vertPathname = NULL;
////	int w, h, d, scale;
////	float stepSize, maxLength;
////
////	ReadConfigureFile(fiber_pathname, &confPathname, &vertPathname, w, h, d, scale, stepSize, maxLength);
////
////	CVector3F *vecs = new CVector3F[w*h*d];
////
////	for (int z = 0; z < d; ++z) {
////		FILE *confFP = fopen(confPathname[z], "r");
////
////		for (int i = 0; i < w*h; ++i) {
////			int id, fnum, bnum;
////			CPoint3F dir;
////			fscanf(confFP, "%d\n", &id);
////			fscanf(confFP, "%d %d\t", &fnum, &bnum);
////			fscanf(confFP, "%f %f %f\n", &(dir.m_x), &(dir.m_y), &(dir.m_z));
////			vecs[z*w*h+i] = dir;
////		}
////
////		fclose(confFP);
////
////		printf("z = %d\n", z);
////	}
////
////	FILE *fp = fopen("\\\\nabla.cs.purdue.edu\\visdata\\Ziang\\CoherenceStructure\\Datasets\\HumanConnectome\\792564_processed\\NEW_DIT_Fibers\\DTI_VEC.raw", "wb");
////	fwrite(vecs, sizeof(float), w*h*d*3, fp);
////	fclose(fp);
////	return 0;
////}

#include "../DWMRICore/Scalar.h"
#include "../DWMRICore/Point3I.h"

#include <vector>

using namespace std;

float threshold_1 = 0.001f;
float threshold_2 = 0.001f;

char *src_pathname = NULL;
char *dst_pathname = NULL;

void PrintUsage()
{
	printf("Usage: RegionGrowing -src <src_file> -dst <dst_file>\n");
}

int ParseArguments(const int argc, char *argv[])
{
	for (int i = 1; i < argc; ++i) {
		if (strcmp(argv[i], "-src") == 0) {
			src_pathname = argv[++i];
		} else if (strcmp(argv[i], "-dst") == 0) {
			dst_pathname = argv[++i];
		} else {
			return -1;
		}
	}

	if (src_pathname == NULL || dst_pathname == NULL)
		return -2;

	return 0;
}

bool FindSeed(CScalar *pSRC, CScalar *pDST, CPoint3I &seed)
{
	uint w = pSRC->GetWidth();
	uint h = pSRC->GetHeight();
	uint d = pSRC->GetDepth();

	for (uint z = 0; z < d; ++z) {
		for (uint y = 0; y < h; ++y) {
			for (uint x = 0; x < w; ++x) {
				if (pSRC->GetScalarFast(x, y, z) < threshold_1 && 
					pDST->GetScalarFast(x, y, z) < 0.5f) {
					seed.x = x;
					seed.y = y;
					seed.z = z;
					return true;
				}
			}
		}
	}

	return false;
}

void RegionGrowing(CScalar *pSRC, CScalar *pDST, CPoint3I seed, const int id)
{
	int w = pSRC->GetWidth();
	int h = pSRC->GetHeight();
	int d = pSRC->GetDepth();

	vector<CPoint3I> buffer;
	buffer.push_back(seed);
	pDST->SetScalar(seed.x, seed.y, seed.z, (float)id);

	while(buffer.size() > 0) {
		CPoint3I p = buffer.back();
		buffer.pop_back();

		float s = pSRC->GetScalarFast(p.x, p.y, p.z);
		

		CPoint3I q;

		if (p.x > 0) {
			q.x = p.x - 1; q.y = p.y; q.z = p.z;
			if (pDST->GetScalarFast(q.x, q.y, q.z) < 0.5f &&
				pSRC->GetScalarFast(q.x, q.y, q.z) < threshold_1 &&
				fabs(pSRC->GetScalarFast(q.x, q.y, q.z) - s) < threshold_2) {
				buffer.push_back(q);
				pDST->SetScalar(q.x, q.y, q.z, (float)id);
			}
		}
		
		if (p.y > 0) {
			q.x = p.x; q.y = p.y - 1; q.z = p.z;
			if (pDST->GetScalarFast(q.x, q.y, q.z) < 0.5f && 
				pSRC->GetScalarFast(q.x, q.y, q.z) < threshold_1 &&
				fabs(pSRC->GetScalarFast(q.x, q.y, q.z) - s) < threshold_2) {
				buffer.push_back(q);
				pDST->SetScalar(q.x, q.y, q.z, (float)id);
			}
		}  

		if (p.z > 0) {
			q.x = p.x; q.y = p.y; q.z = p.z - 1;
			if (pDST->GetScalarFast(q.x, q.y, q.z) < 0.5f && 
				pSRC->GetScalarFast(q.x, q.y, q.z) < threshold_1 &&
				fabs(pSRC->GetScalarFast(q.x, q.y, q.z) - s) < threshold_2) {
				buffer.push_back(q);
				pDST->SetScalar(q.x, q.y, q.z, (float)id);
			}
		}

		if (p.x < (w - 1)) {
			q.x = p.x + 1; q.y = p.y; q.z = p.z;
			if (pDST->GetScalarFast(q.x, q.y, q.z) < 0.5f && 
				pSRC->GetScalarFast(q.x, q.y, q.z) < threshold_1 &&
				fabs(pSRC->GetScalarFast(q.x, q.y, q.z) - s) < threshold_2) {
				buffer.push_back(q);
				pDST->SetScalar(q.x, q.y, q.z, (float)id);
			}
		} 

		if (p.y < (h - 1)) {
			q.x = p.x; q.y = p.y + 1; q.z = p.z;
			if (pDST->GetScalarFast(q.x, q.y, q.z) < 0.5f && 
				pSRC->GetScalarFast(q.x, q.y, q.z) < threshold_1 &&
				fabs(pSRC->GetScalarFast(q.x, q.y, q.z) - s) < threshold_2) {
				buffer.push_back(q);
				pDST->SetScalar(q.x, q.y, q.z, (float)id);
			}
		}
		
		if (p.z < (d - 1)) {
			q.x = p.x; q.y = p.y; q.z = p.z + 1;
			if (pDST->GetScalarFast(q.x, q.y, q.z) < 0.5f && 
				pSRC->GetScalarFast(q.x, q.y, q.z) < threshold_1 &&
				fabs(pSRC->GetScalarFast(q.x, q.y, q.z) - s) < threshold_2) {
				buffer.push_back(q);
				pDST->SetScalar(q.x, q.y, q.z, (float)id);
			}
		}
	}

}

int main(int argc, char *argv[])
{
	if (ParseArguments(argc, argv) != 0) {
		PrintUsage();
		return 0;
	}

	CScalar *pSRC = new CScalar();
	pSRC->OpenScalarFile(src_pathname);

	int w = pSRC->GetWidth();
	int h = pSRC->GetHeight();
	int d = pSRC->GetDepth();

	CScalar *pDST = new CScalar();
	pDST->CreateScalar(w, h, d);

	CPoint3I seed;
	int id = 1;
	if (FindSeed(pSRC, pDST, seed)) {
		RegionGrowing(pSRC, pDST, seed, 1);
		id++;
		//printf("%d\n", id);
	}

	pDST->SaveScalarFile(dst_pathname);

	SafeDelete(pSRC);
	SafeDelete(pDST);

	return 0;
}